HydroData Example’s
Here are some examples of the HydroData Package – April 12, 2018.
Examples in this vingette include
- USGS streamflow
- GHCN
- NHD
- NLCD
- SNOTEL
- Weather Underground
- WBD
- NED
- TIGER roads
- Airports
NOT shown (but available in package)
- Daymet
- Waterbodies
- Reservoirs
- CLD
- PRISM
- Koppen
- NWM
- USGS water use
- NID
Find USGS Gages near UCSB (in a 20 mile square area…)
ucsb.usgs = findUSGS(clip_unit = list("ucsb", 20, 20), ids = T, save = T)## AOI defined as the supplied clip unit.
## Now loading loading CONUS USGS data...## All USGS Data loaded: 6 gages in total## 6 USGS gages found within supplied clip unit## Returned list includes: USGS NWIS shapefile, and list of station IDsThe save = TRUE parameter saves the shapefile to the user working directory (with respect to the query) for later use in R or Arc:
Explore USGS stations
explore(ucsb.usgs, save = T)The save = TRUE parameter saves the map to the user working directory for later reference or sharing:
Find GHCN PRCP records near UCSB
ucsb.ppt = findGHCN(clip_unit = list('ucsb', 20, 20), param = "PRCP", ids = T)## AOI defined as the supplied clip unit. Loading global GHCN data...## 15 GHCN stations found in supplied clip unit## Returned list includes: NOAA GHCN shapefile, and list of station IDsExplore GHCN PRCP stations near UCSB
explore(ucsb.ppt)Get Flow Records from the USGS stations
(shameless wrapper of the USGS dataRetrival function as indicated in documentation):
ucsb.flow = getUSGS(ucsb.usgs$ids)
head(ucsb.flow)## agency_cd site_no Date Flow Flow_cd
## 1 USGS 11119745 1983-10-28 0.81 A
## 2 USGS 11119745 1983-10-29 0.60 A
## 3 USGS 11119745 1983-10-30 0.65 A
## 4 USGS 11119745 1983-10-31 0.70 A
## 5 USGS 11119745 1983-11-01 2.20 A
## 6 USGS 11119745 1983-11-02 0.93 AGet PRCP Records from the first three GHCN stations:
ucsb.ppt = getGHCN(ucsb.ppt$ids[1:3], parameters = "PRCP")
head(ucsb.ppt$PRCP)## agency_cd site_no Date prcp
## 1 NOAA US1CASB0001 2009-09-01 NA
## 2 NOAA US1CASB0001 2009-09-02 NA
## 3 NOAA US1CASB0001 2009-09-03 NA
## 4 NOAA US1CASB0001 2009-09-04 NA
## 5 NOAA US1CASB0001 2009-09-05 0
## 6 NOAA US1CASB0001 2009-09-06 0Visualize the Flows
inspect(ucsb.flow, param = "Flow")## Station 11119745 has 1825 daily FLOW observations.## Station 11119940 has 1820 daily FLOW observations.## Station 11120000 has 1824 daily FLOW observations.## Station 11120500 has 1825 daily FLOW observations.## Station 11123000 has 1825 daily FLOW observations.## Station 11123500 has 1825 daily FLOW observations.Visualize Rainfall
inspect(ucsb.ppt$PRCP, param = "PRCP")## Station US1CASB0001 has 1825 daily PRCP observations.## Station US1CASB0004 has 574 daily PRCP observations.## Station US1CASB0009 has 789 daily PRCP observations.NHD Example
Here we will access a non-point feature
ucsb.nhd = findNHD(clip_unit = list('ucsb', 20, 20))## AOI defined as the supplied clip unit. Shapefile determined. Now loading NHD flowline data...## Trying URL ...## All NHD flowlines loaded: 422 in total.## Returned list includes: flowline shapefileVisualize Flowlines and USGS stations together
explore(list(ucsb.usgs, ucsb.nhd))Land cover/Raster Example
lc = getNLCD(clip_unit = list("ucsb", 10, 10), year = 2011)## AOI defined as the supplied clip unit. Shapefile determined. Now loading NLCD data for 2011.## There is 1 NLCD raster in this scene per year. 1 tiles reuquested...## Downloading raster 1 of 1## Finished downloading raster 1 of 1## Raster number 1 Cropped.Explore Land cover and USGS stations
explore(list(lc, ucsb.usgs))## Warning in colors(.): Some values were outside the color scale and will be
## treated as NAAlternative Means of Describing an Area of Interest
By STATE
Find SNOTEL stations in California
ca.snotel = findSnotel(state = 'CA', boundary = T, ids = T)## AOI defined as the boundary of California. Shapefile determined. Now loading loading Snotel data...## 34 snotel stations found within boundary of California## Returned list includes: snotel shapefile, AOI boundary, and list of station IDsexplore(ca.snotel)Visualize the Snow Water Equivilent for first three stations:
snow.data = getSnotel(ids = ca.snotel$ids[1:3])## 12263 Snotel values for station 301 downloaded (1/3)## 25969 Snotel values for station 356 downloaded (2/3)## 31274 Snotel values for station 1051 downloaded (3/3)head(snow.data)## agency_cd site_no Date SWE PRCP.ACC TMAX TMIN TAVG PRCP.INC
## 1 NRCS 301 1984-09-13 NA NA 32 32 32 NA
## 2 NRCS 301 1984-09-14 NA NA NA NA NA NA
## 3 NRCS 301 1984-09-15 NA NA NA NA NA NA
## 4 NRCS 301 1984-09-16 NA NA NA NA NA NA
## 5 NRCS 301 1984-09-17 NA NA 32 32 32 NA
## 6 NRCS 301 1984-09-18 NA NA NA NA NA NAinspect(snow.data, param = 'SWE')## Station 301 has 1825 daily SWE observations.## Station 356 has 1825 daily SWE observations.## Station 1051 has 1824 daily SWE observations.By COUNTY
Find Airports in Tuscaloosa County, Alabama
al.ap = findAirports(state = 'AL', county = 'Tuscaloosa', boundary = T)## AOI defined as the boundary of Tuscaloosa County, Alabama.
## Shapefile determined.
## Loading Aiport Database## 1 airports found within boundary of Tuscaloosa County, Alabama## Returned list includes: airport shapefile, and AOI boundaryexplore(al.ap)Get WeatherUndergorund Data for airport:
al.data = getWU(airport_code = al.ap$airports$Digit4, year = 2015, month = 1:12, type = 'daily')## Year 2015 Month 1 downloaded.## Year 2015 Month 2 downloaded.## Year 2015 Month 3 downloaded.## Year 2015 Month 4 downloaded.## Year 2015 Month 5 downloaded.## Year 2015 Month 6 downloaded.## Year 2015 Month 7 downloaded.## Year 2015 Month 8 downloaded.## Year 2015 Month 9 downloaded.## Year 2015 Month 10 downloaded.## Year 2015 Month 11 downloaded.## Year 2015 Month 12 downloaded.head(al.data)## Date Year Month Day T_max T_avg T_min DP_max DP_avg DP_min H_max
## 1 2015-01-01 2015 1 1 54 43 31 42 35 29 92
## 2 2015-01-02 2015 1 2 55 50 44 53 47 42 100
## 3 2015-01-03 2015 1 3 68 62 55 66 62 54 93
## 4 2015-01-04 2015 1 4 65 54 42 63 46 32 100
## 5 2015-01-05 2015 1 5 49 40 30 31 25 22 85
## 6 2015-01-06 2015 1 6 56 42 27 33 28 25 92
## H_avg H_min SLP_max SLP_avg SLP_min V_max V_avg V_min Wind_max Wind_avg
## 1 69 46 30.42 30.34 30.27 10 9 5 8 3
## 2 93 86 30.34 30.23 30.14 10 3 1 12 5
## 3 90 87 30.12 30.00 29.92 10 5 1 17 8
## 4 77 54 30.51 30.18 29.97 10 10 4 20 11
## 5 60 34 30.68 30.56 30.42 10 10 8 17 8
## 6 64 35 30.47 30.40 30.32 10 10 7 14 4
## Wind_min PPT_tot Events agency_cd site_no
## 1 10 0.30 Rain WeatherUnderground KTCL
## 2 14 1.61 Rain\n\t,\nThunderstorm WeatherUnderground KTCL
## 3 23 2.47 Rain\n\t,\nThunderstorm WeatherUnderground KTCL
## 4 28 0.37 Rain WeatherUnderground KTCL
## 5 23 0.00 WeatherUnderground KTCL
## 6 19 0.00 WeatherUnderground KTCLVisualize average daily wind velocity
inspect(al.data, param = "Wind_avg" )## Station KTCL has 365 daily WIND_AVG observations.Change timesetep to monthly
inspect(al.data, param = "Wind_avg", timestep = 'monthly' )## Station KTCL has 12 monthly WIND_AVG observations.By a user-supplied shapefile (here all Los Angeles communities subset to Santa Monica)
LA = rgdal::readOGR("/Users/mikejohnson/Documents/GitHub/Communities1/Communities.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/mikejohnson/Documents/GitHub/Communities1/Communities.shp", layer: "Communities"
## with 353 features
## It has 13 fields
## Integer64 fields read as strings: X_CENTER Y_CENTERSM = LA[LA$LABEL_CITY == 'Santa Monica',] %>% spTransform(HydroDataProj)
leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(data = SM,fillColor = "transparent", color = "black", weight = 4)Find all HUC 10 units intersecting Santa Monica
watersheds = findWS(clip_unit = SM, level = 10)## Determining intersecting HUC8 units...## There are 1 HUC8 units in this AOI: 18070104## explore(watersheds)Find Elevation data for the Ballona Creek Watershed (HUC10)
sm.elev = getNED(clip_unit = watersheds$huc_10[1,], res = 1)## AOI defined as the supplied clip unit. Shapefile determined. Now loading 1 arc second NED data...## There are 2 NED rasters in this scene.## Downloading raster 1 of 2## Finished downloading raster 1 of 2## Downloading raster 2 of 2## Finished downloading raster 2 of 2## Raster number 1 Cropped.## Raster number 2 Cropped.## Mosaicing raster...## Returned object contains elevation rasterPlot and ‘analyze’ elevation data for the Ballona Creek Watershed (HUC10)
explore(sm.elev)## Warning in colors(.): Some values were outside the color scale and will be
## treated as NAhist(sm.elev$elev, xlab = 'Elevation (m)', main = paste("Elevation Distribution", watersheds$huc_10[1,]$NAME, "Watershed"), col = 'black')
Find the roads of Santa Monica
sm.roads = findRoads(clip_unit = list("Santa Monica", 2, 2))## AOI defined as the supplied clip unit. Loading TIGER roads database...## There is 1 TIGER file in this scene.## Trying URL ...## All TIGER 1 loaded: 137,781 in total.## Returned list includes: TIGER roads shapefileexplore(sm.roads)